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Standard procedures for local crystal-structure optimisation involve numerous energy and force calculations. 
It is common to calculate an energy-volume curve, fitting an equation of state around the equilibrium cell 
volume. This is a computationally intensive process, in particular for low-symmetry crystal structures where 
each isochoric optimisation involves energy minimisation over many degrees of freedom. Such procedures 
can be prohibitive for non-local exchange-correlation functionals or other ‘beyond’ density functional theory 
electronic structure techniques, particularly where analytical gradients are not available. We present a simple 
approach for efficient optimisation of crystal structures based on a known equation of state. The equilibrium 
volume can be predicted from one single-point calculation, and refined with successive calculations if required. 
The approach is validated for PbS, PbTe, ZnS and ZnTe using nine density functionals, and applied to the 
quaternary semiconductor Cu 2 ZnSnS 4 and the magnetic metal-organic framework HKUST-1. 


I. INTRODUCTION 

The standard operating procedure for computational 
investigations in solid-state chemistry is to begin with a 
crystal structure - obtained either from diffraction stud¬ 
ies or through chemical analogy - and to optimise the 
lattice shape, volume and internal positions to minimise 
all forces. It is from this equilibrium crystal structure 
(athermal for the majority of electronic-structure ap¬ 
proaches) that the full range of material response func¬ 
tions (e.g. elastic, dielectric, magnetic) are calculated. 1 

The optimisation of a crystal structure may require 
hundreds of self-consistent field iterations across a series 
of ionic configurations. 2 The most robust approach to 
optimisation is the calculation of an equation of state 
(EoS) for the material, relating the unit cell dimensions 
to energy and pressure. 3 This is based on a series of cal¬ 
culations at differing volumes, where ideally the shape 
and internal positions are optimised at each point. The 
simplest case is a cubic lattice with high internal sym¬ 
metry, e.g. rocksalt, where the only degree of freedom is 
the volume, and computing the EoS reduces to a series of 
single-point calculations. For a triclinic cell, the lengths, 
angles and internal positions in principle all require op¬ 
timisation. While it is sometimes possible to directly 
optimise the cell volume by simultaneously minimising 
the stress tensor of the unit cell, this approach can run 
into artifacts when using plane-wave basis sets (i.e. Pulay 
forces) unless an iterative procedure is employed. 4 

It has become commonplace to use an ‘equilib¬ 
rium’ crystal geometry, determined using one exchange- 
correlation functional within density functional theory, 
for a ‘single-shot’ higher-level calculation performed to 
give a more accurate electronic structure. This method¬ 
ology has been applied to the calculation of properties 
as diverse as workfunctions, electronic bandgaps, optical 
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properties and defect formation energies. 5-12 The implicit 
assumption is that the qualitative behaviour is insensitive 
to small differences in the local structure. The approxi¬ 
mation will fail where the electronic structure (chemical 
bonding) of a system is poorly described at the initial 
level of theory, e.g. the treatment of Mott insulators such 
as NiO within the local-density approximation (LDA). 13 

In this contribution we outline a simple procedure 
for the rapid volume optimisation (RVO) of crystal 
structures. It takes advantage of the similarity in the 
pressure-volume relationship for a given material between 
different theoretical approaches, here being exchange- 
correlation (XC) functionals within density functional 
theory. Where an EoS is known for one functional, the 
equilibrium volume for another functional can be pre¬ 
dicted with reasonable accuracy using a single calcula¬ 
tion, and further refined following an iterative proce¬ 
dure. The approach has particular utility for studies 
assessing material properties using a range of electronic- 
structure methods, and for studies employing methods 
with high computational cost (e.g. hybrid, meta-hybrid 
and double-hybrid treatments of electron exchange and 
correlation). We validate the approach for four Zn and 
Pb chalcogenides, and demonstrate its utility in describ¬ 
ing the electronic and magnetic structure of one com¬ 
plex semiconductor (Cu 2 ZnSnS 4 ) and one metal-organic 
framework (HKUST-1), respectively. 


II. OUTLINE OF PROCEDURE 

The goal of local crystal-structure optimisation is to 
minimise all degrees of freedom (cell size, shape and posi¬ 
tions) with respect to the total energy of the system. It is 
convenient to employ an EoS based on an energy-volume 
(E-V) curve, where the remaining degrees of freedom 
(i.e. shape and positions) are minimised for each vol¬ 
ume using standard numerical minimisation approaches 
(e.g. the conjugate-gradient method). Kohn-Sham den¬ 
sity functional theory (DFT) 14 is one of the most widely 
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used electronic structure techniques for modelling solid- 
state materials. Most DFT codes provide optimisation 
algorithms for this purpose, e.g. the ISIF=4 setting in 
the Vienna Ab initio Simulations Package (VASP) 15 , the 
cell_dofree=‘volume ; setting in Quantum-Espresso 16 
or the CVOLOPT setting in CRYSTAL. 17 

A superficial resemblance is clear between E-V curves 
obtained with different exchange-correlation functionals, 
with similar shapes but different minima (Figure la). 
The extent of the similarity becomes apparent when 
using pressure-volume (P-V) curves (Figure lb), where 
“pressure” refers to the scalar hydrostatic pressure on 
the periodic system. As this pressure P = — the op¬ 
timal geometries ^ = 0 are now those intersecting the 
P = 0 line. We note that while these still differ depending 
on the chosen XC functional, the P-V curves have sim¬ 
ilar curvature, with the same approximate slopes about 
their zero-crossing points. From these we make our key 
assumption: the slope of one P- V curve may be used to 
estimate the crossing point of another. 

For certain beyond-DFT calculation methods, the 
stress tensor is not computed directly. However, where 
the energy is available the hydrostatic pressure may al¬ 
ways be estimated with a finite difference: 


E(V + 5)-E(V) 
(V + S) - V 


( 1 ) 


The procedure, outlined in Figure 2, is: 


1. Form a P-V curve using one description of the in¬ 
teratomic interactions (method A). This can be 
achieved by fitting an EoS to an energy-volume 
curve. If a system-specific set of classical poten¬ 
tials is available, this would be expected to per¬ 
form very well as they are often fitted to the exper¬ 
imental lattice parameters and elastic properties. 
Within DFT, descriptions of electron exchange and 
correlation within the generalised-gradient approx¬ 
imation (GGA) are suitable, 18 given their low cost 
and the availability of analytical gradients for the 
rapid calculation of forces. Comparative studies 
suggest that PBEsol 19 gives especially good esti¬ 
mates for the lattice parameters and elastic prop¬ 
erties of crystals. 20,21 

2. Calculate the slope about P = 0 for method A. 
This will form our linear approximation. 


m = 


d Pa 

dV 


Pa =o 


( 2 ) 


3. Perform a calculation using a second approach 
(method B), e.g. hybrid DFT with the screened 
HSE06 functional 22 , using an estimated initial vol¬ 
ume; this may be the equilibrium volume (Vo) for 
method A. Convert the resulting stress tensor to a 
hydrostatic pressure Pq. (If no analytical stress ten¬ 
sor is available, use a finite difference as in Eqn. 1.) 


TABLE I: Equilibrium properties of PbS from the 
Murnaghan EoS, fitting over a range of functionals: 
lattice parameter a; unit cell volume Vb; volume 
difference ey from experimental value; Murnaghan EoS 
parameters ko and k' 0 . ko is equivalent to Bulk modulus 
at zero pressure. The experimental lattice constant was 
obtained from low-temperature neutron powder 
diffraction data fitted and extrapolated to zero 
temperature by K. S. Knight (2014). 23 


XC functional 

a / A 

Vb / A 3 

€y / % 

ko 

fco 

LDA 

5.84 

199.01 

-3.47 

65.71 

4.42 

PW91 

5.99 

215.21 

4.38 

55.26 

3.98 

PBE 

5.98 

214.18 

3.88 

54.64 

4.00 

PBEsol 

5.88 

203.43 

-1.33 

61.13 

4.25 

TPSS 

5.96 

211.76 

2.71 

57.33 

4.01 

revTPSS 

5.94 

209.05 

1.39 

57.85 

4.00 

PBE+D2 

5.84 

199.19 

-3.39 

59.93 

5.02 

B3LYP 

6.06 

223.02 

8.17 

53.20 

4.07 

HSE06 

5.96 

210.09 

1.90 

59.29 

4.32 

Experiment 23 

5.91 

206.17 

— 


— 


4. Estimate the corrected volume for method B: 

Vi = Vo+ — (3) 

m 


5. Generate a unit cell with volume V\ (e.g. by in¬ 
terpolating between the previous calculations with 
method A), and recalculate the desired properties 
with method B. 

6. Iterate steps 4-5 until P is acceptably low: 

Vn+% V n • —. P n = f(V n ) (4) 

m 


III. ERROR ESTIMATION 

A. Accuracy of linear approximation 

In this approach, a linear fit is used for the pressure- 
volume relationship: 


Pest — oV + b (5) 

dP e st 


This is by no means a conventional equation of state, 
but may provide a useful approximation when close to 
the minimum volume. The standard definition of the 
bulk modulus, 
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FIG. 1: Energy-volume and pressure-volume curves 
computed for PbS using a variety of DFT 
exchange-correlation functionals. Markers indicate 
calculated values, while smooth lines are fits to the 
Murnaghan equation of state. 


yields the static bulk modulus Bq when evaluated at the 
equilibrium volume Vo- 


B - V dP 
Bo - -V 0 • — 


V = V 0 


( 8 ) 


As we have assumed this derivative to be constant, we 
combine with Eqn. (6) to give a physically meaningful 
expression of our assumption: 


digest _ _ ^>0 

dV ~ a ~ ~Vo' 


(9) 


It is now straightforward to compare this approximate 
EoS with a more conventional form for solid materials, 


estimating an associated error (e). The simplest 
a system with constant bulk modulus. 

dP = _Bo 
dV V 
£ = P — Pest 

de _ dP dPest _ -Bo -B 0 

dE “ dE dV ~ 


V 


V 0 




1(V) = f 

Jv 0 


= s » 


%~ >nV 


1 V 


Vo 




At a typical volume deviation of 5% (Table I): 


e(1.05Vo) = B 0 


1.05Vo - E 0 

V 0 


In 


1.05V, 

% 


= P 0 (0.05 - ln(1.05)) 
where the pressure 


P = 


-/ 

JVq 


s 0 /1.05V4 


case is 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 


and hence the fractional error -p is -2.5%. 


Moving to an improved, while still relatively simple, 
EoS, the Murnaghan EoS adds a parameter, effectively 
giving a linear volume dependence to Bq. 24 Taking its 
derivative form (Eqn. 18), we improve our error estimate: 


dP = _ko 
dE V \ V J 

de = dP dP est = k 0 fVo\ K -B 0 
dV dl dE V V V J Vo 


(18) 

(19) 


We can relate B 0 to the Murnaghan parameters fco, k' 0 by 
finding the slope at Vo- 


Bo 

de 

dE 

<V) 


- v dP 


- V 0u 

v=v 0 Vo k ° 


= kn (l 

l Vo V k o +1 





Yt \ 

yk' 0 +i J 


dV 


= k 0 


V 1 
Vo + K 



k 0 


V 

Vo 


k 0 

= k 0 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 
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FIG. 2: Flow chart for the rapid volume optimisation procedure. 


e(V) = k 0 



= k 0 




(25) 

(26) 


value k' 0 = 5. For smaller values of k' 0 (i.e. closer to 
the fixed-bulk-modulus model), the errors are greatly re¬ 
duced. In any case, the linear approximation appears to 
be sufficiently accurate for a stable optimisation process. 


Plotting these error estimates in Figure 3, we find that 
the error in pressure is less than 10% of the static bulk 
modulus for volumes 10% from the optimum, given a 
material that obeys the Murnaghan EoS with a typical 
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FIG. 3: Divergence of the linear approximation from 
more complex equations of state. The error 
e = P est — Peos is given in units of the static bulk 
modulus. 


We note that \n(x) « x — 1 for x close to 1, and hence the 
residual pressure is approximately linear with respect to 
the error in initial volume estimate. The term p 0,£? ^° ,A 

£>0 ,A Vo ,B 

indicates a smaller linear dependence on the similarity 
of the EoS for method A and method B. Moving to the 
Murnaghan EoS, 


Pi = 


Pi 


Pi 


y 

^0 ,B 

ko,B 

y 

^0 ,B 


o ,B 


V 


k n 


ko,B [ ( Vo,B \ 0 


1 


V 0 , B 


Vo,. 


Vi + P B {V1)\ 

V 0 ,b 


- i 


(33) 

(34) 


\ k O,B 


Vi- 


1 k\ 


Vo,A m ^0 ,B 


0, A K qb 


Vi Vo,Ak o,i 


(Yr 


- 1 


- 1 


(35) 


\Vo,b Vo,Bko,A k' B0 l \ Vi 




- 1 


(36) 


B. Dependence on accuracy of EoS 


Returning to the simplistic EoS of Eqn. 10 


Again, the pressure is dominated by the initial position, 
with a smaller contribution from the difference between 
EoS “stiffness” and volume minima. The non-linearity 
of this relationship follows the non-linearity of the true 
EoS through the exponent k' 0 B . We conclude therefore 
that the performance of the first step is equally sensi¬ 
tive to percentage differences in equilibrium volume and 
bulk modulus between method A and method B. Con¬ 
vergence is impacted by the non-linearity of the EoS, but 


II 

I 

(27) 

not by how accurately this non-linearity is reproduced by 
method A. 

o 

cq 

II 

(28) 

IV. METHODS 



A. Electronic Structure Calculations 


we examine the residual pressure Pi following a single 
step of RYO from an initial volume Y, 


Vi 


Pi 

Pi 

Po,B 


v + =v ._ p B (Vi) Vo ’ A 


m A 

Vi — Po,£ In 

'Vo 9 b 


P. 


0 ,A 


Vqr \ Vo A 


Vi 


Pd,a 


Pd,b In 


Yl 


Vi Po,BVo,A ln fVo,B\ 
Vo,B Pd,aVq,B V Vi ) 


(29) 

(30) 

(31) 

(32) 


Studies have been carried out on the binary chalco- 
genides PbS, PbTe, ZnS, and ZnTe as well as the quater¬ 
nary semiconductor Cu 2 ZnSnS 4 and an organic-inorganic 
hybrid material HKUST-1. 

All DFT calculations on the binary chalcogenides were 
carried out with VASP 25 using the two-atom primitive 
face-centred cubic unit cells. We employed projector- 
augmented wave (PAW) frozen-core potentials 26,27 treat¬ 
ing the outermost s and p electrons of S, Te, and Pb and 
the outermost s, p, and d electrons of Zn explicitly as 
valence. For consistency, we used the LDA PAW poten¬ 
tial set. The PAW potentials are highly transferable and 
tests showed a very weak dependence of the resulting op¬ 
timised electronic and crystal structures. A plane-wave 
kinetic-energy cutoff of 550 eV was employed in all these 
calculations, and the Brillouin zone was sampled with an 
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8x8x8 T-centered Monkhorst-Pack mesh. 28 During elec¬ 
tronic minimisation, the wavefunctions were optimised to 
an energy tolerance of 10“ 8 eV. These parameters were 
found to be sufficient to converge the absolute total en¬ 
ergies to within 1 meV atom -1 , and the stress tensors to 
well within 1 kbar (0.1 GPa). 

The simplicity of the binary systems allowed us to test 
a wide range of functionals, spanning different “levels” of 
approximations to the exchange-correlation potential. 29 
As a baseline, we took the local-density approximation 
(LDA) with the Perdew-Zunger parameterisation of the 
correlation energy. 30 Calculations using the generalised- 
gradient (GGA) approximation were performed with 
the Perdew-Wang 91 (PW91) 31-33 and Perdew-Burke- 
Enzerhof (PBE) 34 functionals, plus the variant of PBE 
revised for solids (PBEsol). 19 To complement this set of 
functionals, we also tested the Grimme D2 dispersion cor¬ 
rection to PBE. 35 “Meta-GGA” calculations were carried 
out using the Tao-Perdew-Staroverov-Scuseria (TPSS) 
functional 36 and the subsequent revision of Perdew et 
al. (revTPSS). 37 Finally, we tested two hybrids, viz. the 
popular HSE06 22 and B3LYP 38 functionals. For each 
material and functional, we calculated an E-V curve by 
adjusting the lattice parameter to yield 21 volumes about 
the experimental 300 K lattice parameters reported in 
Refs. 39 and 40 covering a range of approx. ±5%. We 
note that, as a result of the high symmetry of these sys¬ 
tems, the lattice parameter is the only degree of freedom, 
and thus relaxation of the cell shape and internal posi¬ 
tions was not required. 

For Cu 2 ZnSnS 4 (Section VB), energy-volume curves 
were formed from all-electron DFT calculations using 
the FHI-aims code. 41,42 These calculations employed 
numerically-tabulated atom-centered basis functions (the 
‘tight’ defaults were used, which correspond to expected 
convergence of < 10 meV per atom) and evenly-spaced 
k-point grids. Additional hybrid (HSE06) DFT calcu¬ 
lations and primitive-cell optimisations used VASP with 
the PAW-PBE potential set and a 500 eV cutoff for the 
plane wave basis set. All calculations on Cu 2 ZnSnS 4 sam¬ 
pled the Brillouin zone with 6x6x6 T-centered k-point 
grids. 

For the Cu-based metal-organic framework HKUST- 
1, calculations were again performed with the VASP 
code, considering only the point T in reciprocal space 
due to the large size of unit cell. Owing to the presence 
of the open-shell Cu(II) ion (d 9 configuration) all cal¬ 
culations were spin-polarised, and a range of magnetic 
structures were tested as discussed in Section V C. The 
PBEsol and HSEsol functionals were used along with the 
PAW-PBE potential set. Here ‘HSEsol’ refers to a mod¬ 
ification of HSE06, with PBEsol replacing PBE as the 
local exchange-correlation functional. 43 Due to the com¬ 
plexity of the crystal structure, only three energy-volume 
points were included in the EoS and a single iteration of 
RVO was performed to recover the ground-state HSEsol 
structure. A slightly different procedure was followed in 
this case: a quadratic E-V curve was fitted to the three 


PBEsol points. The initial HSEsol calculation was car¬ 
ried out at the fully-optimised PBEsol point, and the E- 
V curve was followed assuming a constant pressure offset 
to estimate the equilibrium volume for HSEsol. (This 
application was the first chronologically, and led to the 
development of the Murnaghan fitting procedure.) 


B. Implementation 

The RVO method was implemented and tested with 
code written in Python 2.7.3, using the standard li¬ 
brary and Numpy/Scipy/Matplotlib. 44-46 (The code is 
freely available; details in Supporting Information.) 
Non-linear fitting to the Murnaghan EoS uses the 
curve_f it routine in the Scipy Python library, which 
accesses Minpack, an open-source Fortran library. 45 This 
implements least-squares fitting with the Levenberg- 
Marquardt algorithm. 47 Initial guesses of 50.0 eV A 
and 5.0 were used for the k' and k' 0 parameters, respec¬ 
tively. 


V. RESULTS 

A. II-VI Binary Chalcogenide Semiconductors 
1. Simulated application across a range of methods 

For PbS there is a significant range in equilibrium lat¬ 
tice parameters for different exchange-correlation func¬ 
tionals, corresponding to a maximum volume difference 
of over 10%, between the LDA and B3LYP calcula¬ 
tions (Figure 1). Values are tabulated in Table I, and 
compared to a recent low-temperature study by K. S. 
Knight. 23 The computed values are similar, but slightly 
different, to the computational results reported by Hum¬ 
mer et al. 48 Direct bandgaps were also calculated at each 
volume expansion, at the special k-point X for the lead 
compounds and at T for the zinc compounds. (Note that 
PbS and PbTe have another, smaller direct bandgap at 
the L point.) It can be seen in Figure 4 that over the 
lattice-parameter expansion and contraction of up to 5%, 
the bandgaps vary by around 1 eV, with the direction 
of change depending on the structure type and chem¬ 
istry. In this case, using an LDA-predicted geometry for 
a ‘single-shot’ B3LYP calculation would lead to a differ¬ 
ence in bandgap of ~ 0.4 eV compared to that at the 
equilibrium geometry for B3LYP. 

An iterative application of the RVO procedure was 
then simulated from the data. The Murnaghan EoS 
(Eqn. 37) in its integrated form (Eqn. 38) was fitted to 
each E-V curve from DFT calculations. This allowed en¬ 
ergy and pressure to be calculated for arbitrary volumes 
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FIG. 4: Volume-dependence of calculated direct 
bandgaps at T (ZnS, ZnTe) and X (PbS, PbTe) with 
the HSE06 and B3LYP hybrid DFT functionals. 

Results as a function of volume (temperature) for PbTe 
have previously been reported in Ref. 49. The 
behaviour is characteristic of a positive and negative 
bandgap deformation potential for the Pb and Zn 
semiconductors, respectively. The volume expansion 
range of 0.86-1.16 corresponds to lattice parameter 
expansions and contractions of 5% in each dimension. 

Markers indicate the calculated values. 


without carrying out additional DFT calculations. 



The quality of these fits was sufficient for this exer¬ 
cise, with RMS fitting errors of < 1 meV. Fitting pa¬ 
rameters and full data are included in ESI. For each 
“test functional” - “reference functional” pair, the mini¬ 
mum volume (corresponding to the fitting parameter Vo) 
of the reference functional was taken as the initial volume 
guess, and an external pressure calculation modelled by 
evaluating the pressure at this volume using the EoS for 
the test functional. This refined pressure and volume was 
then used as the basis for further iterations. The exter¬ 
nal pressure over successive iterations is shown for PbS in 
Figure 5 for each combination of functionals; convergence 
is rapid with the residual pressure dropping almost log¬ 
arithmically with subsequent steps, typically by a factor 
of ^ 10 3 in three steps. 


2. Comparison with a standard optimisation procedure 

In general terms, a direct optimisation with method 
B will take N optr B steps, each requiring an average com¬ 
puting time £#, to converge to the equilibrium volume. 
Constructing an EoS for RVO using method A requires 
TVeoS ,a optimisations, which, as for method B, take 
N opt:A steps of time We note that in most cases 
TVopt ,b will be larger than N opt ,A, since the direct opti¬ 
misation with method B must adjust the internal coor¬ 
dinates, cell shape and volume, while the optimisation 
with method A needs only to optimise the internal coor¬ 
dinates and the shape. Subsequent application of N^yo 
iterations of the algorithm then requires 1 + TVrvo single¬ 
point calculations using method B, each again requiring 
ts time. RVO is expected to be more efficient than a 
direct optimisation with method B if the following in¬ 
equality holds: 

AeoS,A Aopt,A^A + (1 + ^Vrvo) tB < N opt% bXb- (39) 

The cubic systems considered in this section, for which 
the cell volume is the only degree of freedom, represent 
a special case where N opt: A = 1- We assume that energy 
gradients are available with method B, and that the op¬ 
timisation algorithm would converge in three steps, i.e. 
N opt , B = 3. This is reasonable if a good estimate of the 
starting volume is available, such as a room-temperature 
lattice constant. The inequality simplifies to 


AeoS,aV4 + (1 + Wivo) ts < 3 ts\ 


(40) 
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FIG. 5: Residual pressures following six iterations of volume optimisation of PbS using different “reference”/”test” 
combinations of exchange-correlation functionals. The “reference” refers to the functional used for the EoS 
(“method A”), while the groups of bars correspond to the functionals used in simulated single-point calculations 
(“method B”). Note that the pressure is presented on a logarithmic scale. 


it can be seen that RVO will outperform a direct optimi¬ 
sation if a suitable pressure is obtained after one iteration 

while ^ < ivfc- 

As a concrete example, we compared a direct optimi¬ 
sation of PbS with HSE06 to an optimisation with RVO 
using PBEsol and HSE06 as method A and method B, 
respectively. The initial structure for both optimisations 
was the experimentally-measured room temperature vol¬ 
ume, and an eleven-point EoS for RVO was computed 
about this value using PBEsol. The direct optimisation 
used a quasi-Newton algorithm as implemented by VASP 
(with the input tag “IBRION=l”). Both sets of calcu¬ 
lations were performed on a dual-CPU Intel Xeon work¬ 
station with 12 physical cores and 64 Gb RAM, allowing 
the timings to be compared directly. The comparison is 
given in Table II. 

In this test, a single-point calculation with HSE06 was 
on average 150 times more expensive than a calculation 
with PBEsol; this is both due to the higher complexity 
of non-local hybrid functionals compared to semi-local 
GGA methods, and to the different scaling properties 


with the number of k-points used to sample the Brillouin 
zone. With the force convergence criterion of 10 -2 eV 
A -1 for direct optimisation, the pressure was reduced to 
-0.1 kbar from an initial pressure of 3.42 kbar in three 
steps, taking 4.75 hrs on our test hardware. A single 
iteration of RVO yielded a pressure of 0.15 kbar in 4.18 
hrs, while a second iteration yielded 0.03 kbar in a total 
time of 6 hrs. 

It can be seen from the data in Table II that the direct 
optimisation takes on average less time per force calcu¬ 
lation than RVO; the procedure implemented in VASP 
re-uses calculated wavefunctions to speed up the con¬ 
vergence of the second and third steps. In this case, 
we found that one of the conjugate-gradient electronic- 
minimisation cycles during the first single-point calcula¬ 
tion for the RVO algorithm took some 1500s longer than 
both the other steps in this series and the first step of the 
quasi-Newton volume optimisation, despite the latter be¬ 
ing notionally an identical calculation. This artefact con¬ 
tributes significantly to the cost of the single-iteration 
RVO calculation. 
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TABLE II: Comparison of a direct HSE06 volume 
optimisation and one and two iterations of the RVO 
algorithm in determining the equilibrium volume of 
PbS. For each step, the total time for each step is 
recorded alongside the cell volume and pressure after 
the cycle where appropriate. For the direct 
optimisation, the timings of the three steps are printed 
alongside the total for the complete calculation, so the 
latter includes additional operations such as setup time 
and is slightly longer than the sum of the three 
electronic minimisations. 


Algorithm 

Step 

t / s 

v / A 3 

p / kbar 


1 

6669 

52.21 

3.42 

Direct (HSE06) 

2 

3 

5808 

4509 

52.63 

52.50 

-1.45 

-0.10 


Total 

17038 




PBEsol EoS 

47 




HSE06 1 

8234 

52.21 

3.39 


HSE06 2 

6754 

52.49 

0.15 

RVO 

Total (1 iter) 

15035 




HSE06 3 

6701 

52.50 

0.03 


Total (2 iters) 21736 




Nonetheless, even for this relatively simple test case, 
useful savings in computing time could potentially be 
obtained in practice with RVO. Given the poor scaling 
of computational cost with system size when using ad¬ 
vanced electronic-structure methods, we would expect 
more substantial savings for more complex unit cells. 
This would also be true in systems where direct opti¬ 
misation requires the minimisation of a larger number 
of degrees of freedom, leading a larger number of steps 
with method B, which is the subject of the following case 
studies. 


B. Quaternary Sulfide Cu 2 ZnSnS 4 

Cu 2 ZnSnS 4 (CZTS) is an attractive light-absorbing 
material for thin-film photovoltaics, with a direct 
bandgap and consisting of earth-abundant components, 
which has attracted significant experimental 50 54 and 
computational 55 59 research effort. In the search for new 
materials for solar energy conversion, the prediction of 
accurate bandgaps from first-principles is a serious chal¬ 
lenge and CZTS represents a suitable case for probing 
the effect of crystal structure. 

An initial structure for CZTS in the kesterite phase, 
optimised with PBEsol, was obtained during previous 
work. 60 This was reduced from a conventional 16-atom 
body-centered-tetragonal cell with 14 symmetry to the 
corresponding 8-atom primitive cell using Spglib. 61 A set 
of seven structures was obtained for both cells by mul¬ 
tiplying each lattice vector by a scale factor from 0.97- 
1.03 and performing a local optimisation of the atomic 
positions, this forming the “method A” energy-volume 


curves. In addition to this isotropic scaling, an additional 
set of structures were calculated including optimisation 
of the cell shape (i.e. the tetragonal c/a ratio) for each 
volume point. The iterative RVO procedure was followed 
in order to minimise the pressure and obtain a more ac¬ 
curate electronic structure using the HSE06 functional. 
The results are given in Table III; pressure minimisation 
was rapid in all cases, decreasing logarithmically with 
each step. 

We note that the resulting lattice parameters from 
these calculations, especially those using the primitive 
cell, are very close to both the experimental lattice pa¬ 
rameter a=5.427 A and theoretical a=5.448 A reported 
by Paier et al. following a conventional optimisation pro¬ 
cedure with a variant of the HSE functional. 62 We also 
note a bandgap shift of almost 0.1 eV when the E-V curve 
was provided by a non-isotropic set of primitive lattices. 

This case also highlights the importance of internal 
structure optimisation. After two steps of optimisation 
using the E-V curve further calculations were carried 
out, employing the HSE06 exchange-correlation func¬ 
tional, where the internal atomic positions were relaxed 
while fixing the unit cell. As shown in the table, these 
lead to an increase in the absolute pressure, but also 
a considerable improvement in the bandgap estimation 
compared to experimental measurements. Previous elec¬ 
tronic structure studies have show that the bandgap of 
CZTS is highly sensitive to the S positions, which cor¬ 
respond to deviations away from the ideal tetrahedral 
coordination environment. 56 In the ideal kesterite crys¬ 
tal structure, the metal nuclei all occupy high symmetry 
Wyckoff positions (2 a and 2c by Cu, 2d by Zn, and 2b by 
Sn). However, the sulfur anions occupy the lower sym¬ 
metry 8g positions, with x, y and 2 displacement param¬ 
eters. The change of ~ 0.3 eV in the bandgap following 
further optimisation (Table III) emphasises the impor¬ 
tance of internal relaxations for quantitative studies of 
electronic structure. 


C. Metal-organic Framework HKUST-1 

In 1999, Williams and co-workers isolated Cus(btc )2 
(HKUST-1). 63 Since then this material has been widely 
studied in the field metal-organic frameworks (MOFs), 
with possible applications in catalysis, ionic and electrical 
conductivity, photovoltaics, batteries and gas capture. 64 
First-principles calculations of MOF properties have tra¬ 
ditionally posed challenges for computational chemists 
because they combine large unit cells with complex or¬ 
ganic and inorganic components. 

HKUST-1 features an additional layer of complexity: 
it is composed of Cu-Cu “paddlewheel” inorganic regions, 
where each Cu(II) atom is associated with an unpaired 
electron and each paddlewheel is antiferromagnetically 
(AFM) coupled in the ground state configuration. 65,66 
The magnetic interactions are highly sensitive to the 
Cu-Cu separation. Moreover, previous studies 67,68 have 
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TABLE III: Results from application of RVO to Cu 2 ZnSnS 4 , using the HSE06 functional and a PBEsol-derived E-V 
curve. The unit of pressure P is kbar (10 8 Pa) and volumes V are given in A 3 , a is the lattice parameter in A; these 
are calculated as a mean over the a and b vectors (crystal symmetry is not enforced in FHI-AIMS). E g is the 
electronic bandgap in eV taken from the Kohn-Sham eigenvalues at the T-point. The methods in parentheses refer 
to the process by which the E-V curve was generated; isotropic expansion energies with FHI-aims and 
volume-conserving relaxations with VASP. Iteration “2*” is the data from a final set of calculations, where the 
internal atomic positions are relaxed while maintaining the unit cell from iteration 2. 


Conventional cell Primitive cell Primitive cell 

Iteration (Isotropic expansion) (Isotropic expansion) (Constrained relaxation) 



P 

V 

a 

E g 

P 

V 

a 

E g 

P 

V 

a 

E g 

0 

22.82 

309.12 

5.38 

1.26 

17.49 

155.43 

5.38 

1.23 

17.49 

155.43 

5.38 

1.23 

1 

-1.23 

318.03 

5.40 

1.18 

-0.64 

158.87 

5.42 

1.15 

-1.35 

159.00 

5.44 

1.13 

2 

0.00 

317.56 

5.40 

1.19 

-0.01 

158.74 

5.42 

1.15 

0.02 

158.73 

5.43 

1.14 

2* 

7.46 

317.56 

5.40 

1.49 

10.17 

158.74 

5.42 

1.48 

10.31 

158.73 

5.43 

1.47 


TABLE IV: Results from volume optimisation of 
HKUST-1. Residual pressure P at each step, energies of 
valence band maximum (VBM) and conduction band 
minimum (CBM) with respect to the vacuum level, and 
the bandgap (E g ). All energies are given in eV and the 
pressures is in kbar (10 8 Pa). 


Iteration 

P 

VBM 

CBM 

E g 

0 

-13.98 

-7.5 

-3.7 

3.8 

1 

-1.09 

-7.0 

-3.5 

3.5 


shown that the ionisation potentials and bandgaps of 
porous materials are sensitive to cell pressure and vol¬ 
ume, similar to some of their inorganic counterparts. 
HKUST-1 represents an extreme case, where deviations 
from the equilibrium Cu-Cu distance result in spin flip¬ 
ping and formation of a ferromagnetic (FM) state, which 
impacts the electronic structure. 

Typically, PBEsol-optimised structures agree with low- 
temperature experimental measurements of MOFs to 
within 1 %. This is the case here, and PBEsol also re¬ 
produces the correct AFM state. However, a single-point 
HSE06 calculation on the PBEsol structure yields an in¬ 
correct FM ground-state, as shown in Fig. 6, with an as¬ 
sociated HSE06 cell pressure of -13.98 kbar (Table IV). 
In order to recover an accurate HSE06 crystal structure 
(and the associated correct magnetic structure), a single 
iteration of RVO was required. Notably, there is not only 
a magnetic difference, but also a significant difference in 
predicted electronic bandgap and workfunction. 


VI. CONCLUSIONS 

The rapid volume optimisation (RVO) approach pre¬ 
sented here uses information from an inexpensive energy- 
volume curve to obtain a useful estimate of the optimal 
unit cell volume for a different level of theory. Our focus 


was in bridging between different exchange-correlation 
functionals within density functional theory, but a mea¬ 
sured bulk modulus or classical interatomic potential 
could also be used to construct the reference energy- 
volume data. In sensitive systems the volume change can 
lead to qualitative differences in the electronic and mag¬ 
netic properties. The results depend on the initial volume 
estimate and are relatively insensitive to the accuracy of 
the E-V curve. The RVO method is expected to be com¬ 
petitive with conventional optimisation approaches for 
simple symmetric unit cells, as demonstrated for rocksalt 
structured PbS. For materials such as Cu 2 ZnSnS 4 that 
are sensitive to the atomic positions within the unit cell, 
RVO may form part of the optimisation approach but 
direct optimisation of internal positions is still needed. 
More significantly, it allows for the improved estima¬ 
tion of properties for large unit cells as demonstrated 
for HKUST-1, where conventional optimisation methods 
may be infeasible. As the only property used is the hy¬ 
drostatic pressure, it is possible to employ calculation 
methods which return a total energy without analytical 
gradients by evaluating the energy of a single finite differ¬ 
ence. In this case, an improvement over the “single-shot” 
may be obtained with two additional high-level calcula¬ 
tions and an inexpensive E-V curve; a fourth high-level 
calculation would give an estimate of the convergence. 
We expect that in many cases this will prove an eco¬ 
nomical approach for the application of state-of-the-art 
electronic structure calculations in the solid state. 


VII. SUPPLEMENTARY DATA 

Python code providing a reference implementation 
of this method is available at http://github.com/ 
wmd-bath/rvo; a snapshot as of this publication is 
available at DOI: 10.5281/zenodo.31940. This includes 
a program to generate the plots in this publication 
from the binary chalcogenide data. Calculation data 
has been deposited online with Figshare at the DOI: 
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FIG. 6: (a) HKUST-1 features a periodic array of 32 Cu(II)-Cu(II) paddlewheels per crystallographic unit cell, (b) 
The favoured magnetic structure depends on the Cu-Cu separation: antiferromagnetic (AFM) and ferromagnetic 
(FM) states are accessible, (c) The valence band energy (ionisation potential) is sensitive to the magnetic structure 
(calculated using the procedure outlined in Ref. 67). A ‘single-shot’ HSE06 calculation on the PBEsol structure 
favours the FM state (blue), whilst the corrected structure favours the experimentally observed AFM state (red). 


10.6084/m9.figshare. 1468388. Raw output files from the 
hybrid DFT calculations are made available for CZTS 
and HKUST-1, and energy-volume curves are available 
for all the systems. The full set of graphs and fitting pa¬ 
rameters for PbS, PbTe, ZnS, and ZnTe are also included 
as supplementary data with this paper. 
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